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ABSTRACT 

We discuss the numerical instability arising from the coupling be- 
tween hydrostatic equilibrium and thermal processes in a star. Two 
alternative physical pictures are possible: the heat wave does or does 
not propagate through the adjacent shells in the star in a given time step 
(slow or rapid evolution). Correspondingly, we have two alternative ap- 
proaches to the mathematical formulation. If the physical picture is 
wrong, we encounter a numerical instability. 

In practice, different physical pictures are necessary for a model 
because of a great difference in the time scale of heat conduction between 
the core and the envelope of the star when we compute an advanced phase 
of evolution. After analyzing the nature of the instability, we show that a 
single mathematical scheme is possible which always meets the neces- 
sary physical picture required by rapid or slow evolution in the stellar 
core or envelope. 


I. INTRODUCTION AND SUMMARY 


An elegant, widely-used method fox' the automatic computation of stellar 
evolution was developed by Henyey, Forbes, and Gould (1964), herein referred 
to as the Henyey method. However, the numerical stability of stellar evolution 
computations has not yet been analyzed sufficiently. The problem of the com- 
putation of stellar evolution is characterized as a mixed initial-boundary value 
problem with four simultaneous differential equations , two describing hydro- 
static equilibrium and the other two the thermal process. In the present paper 
we shall discuss the most serious numerical instability which arises from 
coupling between hydrostatic equilibrium and the thermal process. 

a) Explicit and Implicit Methods 

We shall consider a method of obtaining the stellar structure at time t , 
assuming that the preceding structure at time £ - A t is known. The star is 
divided into a suitable number of shells. There are two alternative ways of 
solving the problem. One is called the explicit method, in which any quantity 
at t is expressed explicitly in terms of the quantities at t. - At.. The other is 
called the implicit method, in which the quantities at 1 are expressed implicitly 
with respect to the quantities at jt - A jt, and some elimination method is required 
to obtain the quantities at .t . The Henyey method is an implicit method in this 
sense. 

There are mixed expressions; for example, when the thermal process is 
expressed explicitly but the hydrostatic equilibrium is expressed implicitly. 
Examples of such expressions were given by Schwarzschild and Selberg (1962), 


by Rakavy, Shaviv, and Zinamon (1967), and by Murai, Sugimoto, Hoshi, and 
Hayashi (1968). In these examples, the entropy density, s (M r , t) t is calculated 
by using quantities determined only at t - At.*, where M r denotes the mass 
inside the radius r.. The notation in the present paper is the same as that given 
in the textbook by Schwarz schild (1958) unless otherwise indicated. By using 
s (M r , t), we can solve the two differential equations of hydrostatic equilibrium 
for time t separately from the equations for the thermal process. The heat 
flux, L r (t), is calculated from the temperature distribution after the structure 
at t, has been obtained. In the following discussion, it is an essential point as 
to whether the thermal process is expressed explicitly or implicitly; the expres- 
sion for the hydrostatic equilibrium does not matter. Thus we shall call a method 
"explicit" whenever the thermal process is expressed explicitly. 

b) Stability Condition for a Purely Thermal Process 

The stability of a heat-flow problem in a solid body has been well studied 
(see, e.g., Richtmyer and Morton 1967). The stability condition of the explicit 
method can be written as 


2At ^ r h (Ar) [explicit] , 


( 1 ) 


Rakavy et aj_, (1967) hitroduced a damping term in quasi-time, t q . Thus the entropy density 

was expressed as s (M f , t; t ), with s (M f , t) » lim s (M r , t; t ). They calculated s (M r , t;t ) 

q t -.w s s 

q 

by using s (M f , t; t q - At q ) and s (M f , t - At). This was an iteration of the explicit method 
for the thermal process, since the mechanical equilibrium was solved separately by assuming 
the entropy denisty. 
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where the time scale of heat conduction through a shell of radial width A r is 
given by 


T h (Ar) 


’rod 


Bs 

r) In T 


K 


f) (A r)' g 


( 2 ) 


In the above equation, s rnd denotes the entropy density of radiation (4aT 3 /p), 
On the other hand, the explicit method is stable for any time step in a purely 
thermal problem. 


c) Stability Condition for a Coupled Hydrostatic-Thermal 

System (Real Star) 

For the explicit method the stability condition is the same as in relation 
(1), since the hydrostatic equilibrium is solved separately from the thermal 
process, as discussed in § I, a above. For the implicit method, on the other 
hand, the condition that the system of four differential equations be stable is 
expressed roughly by 


At £ r h (r) [implicit] , (3) 

where r h (r) denotes another time scale of heat conduction expressed by equa- 
tion (2) with A r replaced by jr . We call a system of differential equations un- 
stable when it has a branch of a solution growing sc strongly that it is practically 
impossible to compute all of the independent solutions numerically (Wendroff 
1966). This point will be discussed in more detail in § II. 


cl) Problem Encountered in Practice 


We shall now consider the computation of stellar evolution through a phase 
with a reasonable number of time steps. As may be seen from the stability 
conditions (1) and (3), the explicit method is suitable for rapid evolution (At small 
compared withr h (Ar)), c.g. , for evolutionary phases with extensive neutrino 
loss or thermal instability, but, on the other hand, the implicit method is suit- 
able for relatively slow evolution, e.g., in phases of nuclear burning with neg- 
ligible neutrino loss. 

The time scale of heat conduction varies greatly through the star, especially 
in a red giant star with a helium zone and hydrogen envelope.* For a reasonable 
division into shells, r h (A r) has its smallest value at the outermost shell in radi- 
ative equilibrium, even where condition (1) must be satisfied in the explicit 
method. (In a convective region, r h can be considered to vanish.) For a reason- 
able value of At, in some cases, we encounter a contradictory requirement that 
the implicit method is necessary in the outer region but the explicit method is 
required in the central region. Typical examples are the phases of and near 
carbon burning. 


e) Stability of Difference Equations and Conclusion 


Fortunately, such a contradictory requirement can be satisfied, as will be 
shown in § HI. We shall look more carefully into the behavior of the system of 


For example 7^ (A r) for a unit scale height of pressure, h_, is 2 x 10^ sec at the central region 
(r =* h), 8 x 1011 sec at the helium-burning shell, and 3 x 10 10 sec at the hydrogen-burning shell 
of a 15 M 0 star just before the carbon-burning phase (stage 5, Table 7-6 of Hayashi et ah (1962)). 



difference equations of the Henyey method, In the limit of an infinitesimal time 
step but keeping spatial step finite, the difference equations no longer approxi- 
mate the original system of four differential equations.* It is shown that, with 
a properly formulated Henyey method, wo can obtain physically significant 
branches of solutions by avoiding the unstable branches of solution which the 
original system of four differential equations has in the above limit. Solution 
by such a method can be expressed by an expansion in which the leading term 
is a solution of the explicit method and the first-order term represents the 
coupling between a change in the heat flux and a change in the hydrostatic equi- 
librium. Since this method is essentially a form of the Henyey method, it is 
stable for slow evolution. Thus, we can solve stably both the core and envelope, 
as well as both rapid and slow evolution, in a single scheme of computation. 
Application of this method to a physical problem will be given in a separate 
paper (Sugimoto 1969). 


H. NATURE OF THE DIFFERENTIAL EQUATIONS 
FOR STELLAR EVOLUTION 

a) Differential Equations for Stellar Structure 

The structure of a star in quasi-static equilibrium is determined by the 
following four differential equations: 

B In p _ ^ 

3 In M r 477T 4 p 

‘k 

In order to approximate the original differential equations, both the time step and the spatial 
step must become infinitesimal simultaneously, keeping At/Ar finite. 
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(4b) 


3 in T 1 i In p 
3 In M f n + l 3 In M r 


B In r 

0 In M r 47 /t^ p 


(4c) 


where 


BL„ 


V) I n M„ 



n + 1 


min 


1 1 
(n + ^)ad + ^r»d 


(4d) 


(5a) 


B 1 — - 2. — (5b) 

(n + l )md 16 17 0CG T 4 M r 

and l/(n + l) d is the adiabatic temperature gradient. For the sake of definite- 
ness in the following discussion, we take the dependent variables as In p (pressure), 
In T (temperature), In r, and L r ; these will be denoted as y 1( with i - 1, 2, 3, 
and 4, respectively. Hereafter we shall denote a matrix by a capital letter, a 
vector by a lower case letter, and their elements with subscripts. The inde- 
pendent variable, In M r , will be denoted by x . The nuclear energy generation 
rate, e n , the energy loss rate by neutrinos , e v , and the opacity, k , are usually 
expressed in terms of p andT . The density, p, the entropy, and the adiabatic 
temperature gradient are in principle expressible as functions of j> and T. 
Equations (4a) -(4d) are rewritten as 

|2-=?(y.x)- (6) 

O X 
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Of course, transformations of the dependent variables can be incorporated, but 
an essential point in the following discussion is that the equations for hydrostatic 
equilibrium involving y l and <p 3 do not contain y 4 , There are two boundary 
conditions at the center: r * L r * 0 at M r 0, There are twv other boundary 
conditions at the surface: p * T * 0 at M p * M (total mass of the star), when the 
surface region is in radiative equilibrium, 



b) Relaxation Metho d 

The relaxation method was first introduced by Henyey et al. (1959), for 
solving problems in stellar evolution. The time derivative in equation (4) Bs/Bt , 
is replaced by the forward difference, ie M by { s (t ) - s (t - At)}/At. Except 
for this term, equation (4) does not contain any time derivative. We shall hence- 
forth write d/dx instead of B/Bx . We assume a trial solution y<°> and assume 
that 

y B y(°> f §y (7) 

satisfies equation (6). Assuming that £y is a small correction to y(°> , then 
substituting equation (7) into equation and linearizing it, we obtain 

4^«A$y + fa, (8a) 


* 



(8b) 


where A i(j - (B<p./By.) (0) and b. = <p. (y <°>, x) - dy i <°Vdx are numerically 
known functions of x . After having solved equation (8a) and obtained y by (7), 
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y< 0) is replaced by y . The whole procedure is repeated until by becomes suit- 
ably small. Thus the problem is how to solve equation (8). 

c) Nature of the Differential Equations 
for the Correction Term 

We consider only a region in radiative equilibrium, since in a convective 
region all of the A. _ 4 vanish and y 4 is separable from the others. If a finite 
efficiency of the convective heat-transport is taken into account by means of 
mixing-length theory, for example, it can be treated in the same way as the 
radiative equilibrium in the sense that l/(n + 1) depends upon L f . Important 
quantities in the following discussion are written as 



We consider the computation of a phase of relatively rapid evolution with N 
time-steps. Where the heat conduction is negligible, the left hand side of equa- 
tion (4d) should be small compared with each term in its right hand side. There 
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is convection wherever the nuclear energy generation is large. Such a convective 
region does not contribute to the numerical instability. Therefore the neutrino 
loss, & v . approximately equal to {s (t) - s (t - At) } T/At. The change of en- 
tropy density, |s(t) - s(t - At) | , is at most about (3s/3 In T) p /N, and less than 
this for a shell in which the neutrino loss is relatively small as compared with 
other shells in the star. Thus the terms proportional to (3s / 'd in p) T and 
(3s /3 InT) are most dominant in A 41 and A 4 2 , Taking into account only these 
terms, we have 


n =(Ai,4-V») ! 



T l/ r > 

At 


U 


( 10 ) 


where Usd In M r /d in r is three times the ratio of the density to the mean in- 
terior density as seen in equation (4c). In the limit of infinitely rapid evolution 
or an infinitesimal time step, fi is infinitely large, while A /A. „ remains 
finite. 

We consider a range of x , where A can be considered to be practically 
constant. The nature of the solution for the homogeneous part of equation (8a) 
is understood by the secular equation, 

| A - K | =0. (11) 

Taking into account equation (8b) , it is easily shown for large H that two eigen- 
values are given by 

a 4 = - s a + o <n°> . (12) 
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The other two are given by 


\ a \ k < 0) + 0 (Cl" 2 ), k B 1, 2, 


( 13 ) 


where \ k (e) denotes eigenvalues of the equations 


d 

' dx 



A i,i 


i ^ 1| 3 1 



(14a) 


(14b) 


These are the eigenvalues of the two differential equations for hydrostatic equi- 
librium, i.e., equation (14a), together with a given distribution of entropy as a 
subsidiary condition (14b); they are the eigenvalues of the explicit method. Only 
these two eigenvalues are physically significant in the limit of large Cl. 

The general solution of equation (8a) contains strongly growing branches, 
exp (jfcHx ), and the differential equations themselves are unstable. Thus, it is 
impossible to obtain independent solutions by numerical integration. Even if one 
has started from different sets of initial values, the unstable branch overcomes 
the other branches after some steps of integration, and these solutions are no 
longer practically independent of each other (Wendroff 1966).* 

A measure of the criterion for equation (8) to be practically stable is given 
by using equation (10) in the form 

* ^surface 

I = n dx ~ (r h (r)/A t } 1/2 d in r ~ 0 (1) . (15) 

J Jr B 0 

Even when we start with properly selected initial values so as not to contain an unstable 
branch, approximation with a finite step of integration introduces the unstable branch. 
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The stability condition (3) given in § I, c is another expression of the above con- 
dition. This means that a heat wave must propagate throughout the star in a time 
At, as required from the physical picture of equation (4). It is to be noticed that 
the integral in the above equation applies only in the radiative region. 


in. DIFFERENCE EQUATIONS 


We consider only the homogeneous part of equation (8) with a difference 
equation of the following type: 


- 8 y 1 ^ 1 + 8 y*5 + ^ yS. A^ ( . Ax k 8 y 1 ? 


j*i 


4 

+ Y -1 (1-/3,) AH 1 Ax k 8y^ + 1 * 0, ( 16 ) 

j = i 

where the superscript denotes the spatial mesh points, k = 1, 2, . .. ,K-1 and 


Ax k = x k i- 1 - x k . 


(17) 


We are concerned only with the coupling between hydrostatic equilibrium and 
thermal process so that we put /3 t = JS 3 = - . 


a) Henyey--Type Eli mination 
We introduce new independent variables, 

s v \ = s y \, i = 1, 2; S-n\ = S + 1 , i =3,4. ( 18 ) 
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aaa assasaa 


Equation (16) then becomes 


P k 8 7 ? k “ 1 + Q k 8 r) k + R k 8 7? k + 1 « 0, (19) 

Of course the first two columns of the matrix P k and the last two columns of R k 
are zeros, introducing another matrix, P k , of which the last two columns are 
zeros, and which satisfies 

8 a P k g 77 k , (20) 

the solution of equation (19) is given by 

r k + 1 s (~p k r k - Q k ) _1 R k , (2ia) 

& r) k = r k + l r k + 2 — r K st] k . (2ib) 

We have two independent choices of 8 r) K ., since the choice of 8 773 and 8 77$ is mean- 
ingless. On the other hand, f! . (i = 1,2, j = 1,2) does not enter into equation 
(21a). Thus we have two independent choices of T l , for which the determinant of 
the submatrix consisting of r! . (i = 3,4, j = 1,2) does not vanish. Consequently, 
v/e have two independent sets of F k . Thus we have four particular solutions.* 

b) limit to a Small Time Step with Fixed Mesh-Points 

We now consider the case with a large value of Q . If we solve equation (16) 
with mesh-points satisfying co = flAx «1, the difference equation (16) approximates 

Usually, the two boundary conditions at the center are incorporated in equation (19) for k = 1, so 
that P vanishes and the choice of is meaningless. Then, we obtain only one set of T k , i.e., 
two particular solutions which satisfy the two boundary conditions at the center of the star. 
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the differential equation (8) and we pick up the unstable branch. Thus, we ex- 
amine equation (21) for the case of <v » 1, i.e., for rapid evolution with a finite 
number of mesh-points. The result depends upon the detailed formulation of the 
difference equations as well as upon the method of elimination.* 7t is difficult to 
show generally what the solution (21) approaches in the limit of a large value of o> . 
However, numerical experiments show that it approaches the solution of the ex- 
plicit method when we choose 


or 


/? 2 

fit 


0 and /3 4 s 1 , 

(22a) 

1 and j 3 4 = 0. 

(22b) 


We shall examine mathematically the first case, only because a mathemati- 
cal proof is easiest in this case. We drop the superscript when no confusion is 
anticipated. Taking account of equations (8b), and (22a), the matrix elements 
are written as 

w S -P k r k - Q k 


For example, if we assume two sets of values for Sy 1 , which satisfy the boundary condition at 
the center of the star, we can solve equation (16) numerically for the two particular solutions. 
Numerical experiments show that these two solutions are not practically independent at the 
stellar surface, reflecting the nature of the differential equation (8). 
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“ P l,3 P 3,l “$1,1 
“$2,1 > 

" P 3,3 ^3,1 “ ^3 , 1 


P l,3 * 3,2 "Ql ,2 ’ ~Ql,; 


~ ^2 . 1 ' 


“ P 3 , 3 * 3,2 " ^3 , 2 


"Q2.3* ~Qq 


-Qa.a. 0 


\- P 4,3 r 3, 1 “ r 4,l -04.1* “ P 4,3 * 3,2 ~ **4.2 ~ $4,2 » 1 J 

We split the matrix W into two parts so that | W |= | W (c) | + | | : The matrices 

and w rf > denote the same matrix as W but the fourth rows are replaced by 

(-Q 4 ,i» -Q 4 , 2 * °* °> and (“ p 4 , 3 r 3 ,i “ r 4 ,i' - p 4 , 3 r 3 , 2 - r 4 , 2 ’ 0, 1 ), respectively. 
Hereafter, we shall not write Ax explicitly, since it is fixed. The element, -Q 2 t , 

appears if we do not substitute equation (4a) into (4b). It should be noticed that 
Q 4 2 Q 2i 4 is equal to w 2 and is invariant for a scale change of S 77 j . 

We assume that Q 2 , 4 r 4 , 1 » Q 2 , 4 P 4 , 2 and the other r 1 j are sma11 compared 
with w 2 as should be the case for k « X. Denoting the co-factor of W £ . } by | W. ( i \ , 
we have 



Bince }w£®^| is a determinant which appears in the explicit method, it does not 
vanish insofar as the explicit method is stable. Thus,Q 2 4 | is of the order 

of co 2 and the quantity in the curly brackets in equation (24) is of the order of co" 2 . 


We denote a matrix as W (R; i, j )which is the same matrix as W but with the 
i-th column of W replaced by the j-th column of R. Using equations (21a) and (24), 
and noting R 41 = R 4 2 = 0 because of equation (22a), the other matrix elemei^s of r 
can be expanded in w" 2 as 
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pk + l . 
i. J 



. | f!w rf) CR;i.i) a , 4 l 
[!§(«» (R;i,j) 2 ,J 





where the quantities in the curly brackets are of the order of of 2 , Thus, all of 
r i,V and Q 2 ,V r 4 , + j 1 remain finite since is finite. 


c) Physical Interpretation 

The leading terms of equation (25) and of Q 24 r$ + 1 inequation (26) contain 
neither the (2,4)-element, T 4 nor any co-factor of the (4 ,4) -element. More- 
over, the leading term of equation (25) does not contain the (2,j) elements. Thus, 
the leading terms are solutions to be obtained by the explicit method, i.e., by 
letting P 4 4 and Q 4 4 vanish. In spite of a large value of a> , these are correct 
solutions of equation (14a) with the subsidiary condition (14b), since 
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A. ( j Ax (i,j = 1,3) can be small compared with unity. The by A is calculated from 
d&y 2 'dx, Thus, the leading terms describe solutions for given distribution of 
entropy density as discussed in § I, a. The first-order terms in m~ 2 describe effect 
of heat flow, and thus these are proportional to the co-factors of the (2,4) or the 
(4,4)-elements. 

Thus, the mathematical structure of the above scheme corresponds to a good 
physical picture that the heat wave does not propagate in a given time interval 
through a shell in the limit of large co. On the other hand, when co is small, this 
scheme reduces to the usual implicit method and thus corresponds to a good 
physical picture that the heat wave propagates well throughout the star. The case 
of an intermediate value of co will be discussed in § III, d. 


We must now discuss the number of independent solutions and the boundary 
condition for a large value of co throughout the star. We easily find that 


pk pk 

1 1,1 1 1,2 


pk pk 

2.1 1 2,2 


= o <y 2 ), 


(27) 


which means that the first two columns of p k p k + 1 , . , p K are not independent in 
the limit of infinite co. Thus we have only two degrees of freedom in the choice 
of a solution, i.e., only in the choice ofr 1 . This corresponds to the fact that we 
have only two differential equations of the explicit method in this limit. One of 
the boundary conditions in this limit at the center, L r = 0, has nothing to do with 
the hydrostatic equilibrium. One of the boundary conditions at the surface be- 
comes a relation for the entropy at the outermost point. Thus we have only two 
boundary conditions — one at the center and another at the outermost point - 
which are satisfied by using the two independent solutions. If co may be 
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regarded as Infinite in the shells for k» 1,2, , , , , k 0 but small in the shells 
for k = k 0 + 2, . . . , K, we have four independent solutions in the outer shells 
but these degenerate into two Independent solutions in the inner shells. Thus 
the two boundary conditions at the outermost shell and one boundary condition at 
the center can be satisfied correctly. Although we have not discussed the in-* 
homogeneous term in equation (8), It does not alter the above discussion, 

d) Discussion and Comments 

It is difficult to treat the case when w takes an intermediate value. How- 
ever, the following discussion will give some idea regarding a practical case. 

We must consider a region where 1 < 10, In this region the strongly- 

growing branches (equation [12]) of the differential equation cannot be solved 
exactly, since w = HAxis too large. Thus the integral (15) must be replaced 
by aP , where n is the number of such shells, If a> n is smaller than 10 H , we 
can extract the physically significant branches (equation [13] ) by using the double 
precision computation. Although branch (12) is not solved exactly, there will be 
no problem, since dy 4 /dx will be small enough compared with | A 4 ^ | in such 

a region. 

We have not discussed other choices for equations (18) and (22). For ex- 
ample, the choice of (S 77 k) as (Syj* , Sy 4 , SyJ + 1 , 8y£ + 1 ) prevents subtraction in 
calculating SyJ from other variables by using equation (20). However, it is to 
be noticed that even in the choice of equation (18), the rounding-off error does 
not accumulate. For a phase or a region in which is small enough, fi 2 = /5 4 = i- 
or other weighting gives a better approximation and, a more rapid convergence 
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of the relaxation computation than equation (22), However, it Is only a matter 
of programming technique to take this into account, 

In summary, the Henyey method is a good mathematical scheme based on 
good physical pictures of both slow and rapid evolution if it is properly formu- 
lated, We can solve stably the radiation flux, L r , as well as other variables, 
even for rapid evolution. The stability of L r is a good measure of whether the 
method has beon formulated in a given case so as to represent a good physical 
picture. 
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